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Background 


The purpose of this Grant is to investigate the possibilities of using hp-time finite elements to 
solve problems of rotorcraft trim and periodic response. Although conventional time-marching and Fourier 
methods have been somewhat successful, they do not always converge. Thus, methodologies need to be 
studied. The class of hp methods in time has the advantage that it is applicable even to unstable systems 
and easily converges on discontinuities in the solution. 

Personnel 


During the duration of this Grant, the Principal Investigator (Dr. David Peters) spend 0.9 man- 
months of effort; and Graduate Research Assistants Ay Su, Lin-Jin Hou, and Yi-Ren Wang spent 2.57, 
1.00, and 0.75 man-months, respectively. The work of Mr. Su and Mr. Wang was primarily in providing 
flap-lag response for comparisons. Mr. Hou, however, did all of the finite-element coding. 

Results 


The results of this work include trimmed and untrimmed cases, flap and flap-lag cases, and are 
for both displacement and mixed formulations. The attached paper, to be presented at the AHS 
International Specialists’ meeting on Rotorcraft Basic Research, summarizes the results. In addition, an 
oral presentation was given at the DAMUIBS workshop in September, 1990, at NASA Langley. 
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ABSTRACT 

We have studied finite elements in time as an alternative strategy for 
rotorcraft. trim problems. The research treats linear flap and linearized flap- 
lag response both for quasi-trim and trim cases. The connection between 
Fourier series analysis and hp-finite elements for periodic problem is also 
examined. It is proved that Fourier series is a special case of space- time 
finite elements in which one element is used with a strong displacement for- 
mulation. Comparisons are made with respect to accuracy among Fourier 
analysis, displacement methods, and mixed methods over a variety param- 
eters. The hp trade-off is studied for the periodic trim problem to provide 
an optimum step size and order of polynomial for a given error criteria. It is 
found that finite elements in time can outperform Fourier analysis for peri- 
odic problems, and for some given error criteria. The mixed method provides 

better results than does the displacement method. 

1 Graduate Research Assistant, School of Aerospace Engineering. 

2 Professor, School of Aerospace Engineering. 
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INTRODUCTION 


Finite elements in the time domain have recently come into widespread 
use. Versions exist in a time-marching framework, in a Galerkin framework, 
in a Rayleigh-Ritz framework, and in a mixed formulation. Furthermore, 
recent hp work has demonstrated that the convergence of such methods is 
highly enhanced when one uses a weak formulation with test functions chosen 
so as to meet certain precise mathematical criteria on the continuity of the 
bilinear operator [13]. Thus far, these fast-converging hp methods have only 
been used on transient response (Floquent theory) [13] and on boundary- 
value problems (in optimal control) [14]. In this paper, we use this method 
on rotorcraft trim problems. 

For hingeless rotorcraft dynamics, aeroelastic motions of the rotor blades 
and of the fuselage are coupled in forward flight thus makes the response 
problem fairly difficult. Analyses of pitch-flap and pitch-flap-lag dynamics 
in hover and flap dynamics in forward flight have been well documented [3] 
[6] [10], and rotorcraft trim is usually one of the most difficult aspects of the 
solution. In brief, to trim, one must find a periodic solution to a nonlinear, 
periodic-coefficient, differential equation subject to side constraints that cer- 
tain force and moment balance equations must sum to zero. There are certain 
free (or trim) parameters in the problem that must be chosen so as to meet 
these side constraints, but it is not easy to find either the periodic solution 
or the unknown controls. Although conventional time-marching and Fourier 
methods have been somewhat successful, they do not always converge. 

Hamilton’s law of varying action has been applied to time finite elements 
[1] [2] [4]; recently, Borri et al [5] and Peters and Izadpanaph [13] provided 



a powerful alternative by using Hamilton’s weak principle (HWP) to the nu- 
merical solution of ordinary differential equations in the time domain. Here, 
we apply Hamilton’s weak principle for time finite elements to the problem 
of rotorcraft trim and periodic response. The work begins with the question 
of shape functions and compares the use of integrals of Legendre polynomi- 
als as shape functions to the use of Fourier series; since the latter enforce 
momentum periodicity (which is a natural boundary condition) they force 
a strong rather weak formulation. However, we need to firmly establish the 
connections between Fourier analysis and hp-finite elements for the periodic 
case. 

Second, the trim constraints are treated as formal side constraints in the 
variational form of the problem. This turns out to be a benificial by-product 
of the hp methodology - namely that the unknown controls can augment 
the unknown coefficients in a formal way. In addition to the displacement 
method, the mixed method is also implemented in which displacement and 
momentum have seperate expansions. Finally, the inplane degrees of freedom 
is added. In terms of numerical issues, the hp trade-off is investigated; and 
the work herein will then show that how the hp trade-off is optimized for 
both methods, as advantages are verified for accuracy and efficiency. 

• Linear Flap 

For the linear flap problem, the equation of motion is derived for an 
idealized model of a hingeless rotor rigid blade. Aerodynamic forces are 
based on linear quasi-steady strip theory; small angle assumptions are made 
for induced angles and for the flap angle. The induced inflow is assumed 
to be uniform and is obtained from simple momentum theory. Also, the lift 


curve slope o, and blade chord c are assumed constant along the blade. 

The equation of motion for flapping can then be obtained from moment 
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equilibrium 

0 +C(i)/3 + K(i) 0 = F(t) (1) 

where 

C(t) = | (1 + ^ sin i) 

K(t) = p 2 + cost + /X 2 sin 2t) 

O O 

g 

F(t) = -[(0 O + 0, sinf + 0 e cosi)(l + -/xsin t + p 2 - p 2 cos2f) 

8 d 

4 - 4 

— A(- -f 2/i sin t) — 4>(1 + -/xsint)] 

3 v 

i = v = nt, ( ) = ^( ) 

FOURIER SERIES ANALYSIS 


To implement a Fourier analysis, we express the flapping angle (3 in terms 
of a harmonic series with JVj, harmonics and the number of terms N — 2 Nh + l 

0 = £¥>*(<>>■ 

*=i 

= qi + q2 cos t + <73 sin t + <74 cos 2 1 + <75 sin 2 i + ... 


We then apply Galerkin’s scheme, multiply equation (1) by test functions <fj, 
which are also harmonics, and integrate. We then have 


V f (<fi + Cpi + K^i) qWjdt = f Fifj dt j = 1, 2, . . . , N 
fr{Jo Jo 

Integration by parts and use the periodicity of the tpi yields 

N r 2r - f 2 " 

£ / - c Vi<Pj - K <Pi<Pj) qidt = - Fifijdt 

f^ J ° J ° 


( 2 ) 


( 3 ) 



Now, we consider for comparison the alternate approach of using finite 
elements in time by virture of Hamilton’s strong principle 

6 f\T - V)dt + SW ne dt = 0 (4) 

Jo Jo 

where 

T = U\ V = \k&\ 6W ne = (F-C(3)6(3 

for the current problem. Substitution of the above into equation (4) yields 

[ 2ir {i36(3 - K(36(3 + F)3 - C(36(3)dt = 0 . (5) 

Jo 

Let 

0 = y. &$)* , w = n 

i=i j=i 

where the <fj must be periodic, since equation (4) must be true for arbitrary 
We will then obtain equation (3) as deriverd from Fourier series. The 
governing differential equation can also be recovered from equation (5) b} 
integration by parts and use of the periodicity property on momentum. Thus, 
this proves that the Fourier series approach is a special case of finite element s 
in time in which harmonics are used as trial and shape functions for one 
element. 

For the quasi-trim problem, the collective and cyclic pitch angles are 
given. 8 = ~8 9 sin so as to suppress approximatly first harmonic 

cyclic flapping, 3 C , f3 s ^ 0, and thereby minimize rotor hub moments. For the 
trim problem, the trim parameters 8 o, 8 C , 8 S are treated as unknown controls 
to be found such that first harmonic flapping is identically zero. 

The equations for these trim conditions, which can be taken as formal 


side constraints, are 
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— f j3costdt = 0 (7) 

7 r Jo 

— [ f3sinidi — 0 (8) 

7T Jo 

Due to the orthogonality of triogonometric function, equations (7), (8) 
yield for Fourier analysis 

<?2 = 0 , q 3 = 0 

Equation (6) can also be replaced by an alternative form in which thrust is 
specified 

— — f Cr{i) dt = Cto 

2t Jo 

where Cy(f) is obtained by blade element theory 

C T (t) = [ + fi sin tf 6 - (x + n sin t)(\ + xfi + fi/3 cos t)]dx 

J o 2tt R 

For the trim condition, the minimum order necessary for the Fourier series 
approach is first harmonic cyclic flapping (<ji, qz-, 93); since with use of 91 only 
as a one-term approximation (i.e., (3 = flo), equations (7) (8) will be identi- 
cally zero which results in a singular matrix. Physically, three polynomials 
are required for Fourier series ( Nh = 1) because of three trim constraints; and 
the periodicity conditions are automatical^ satisfied. Similar requirement on 
the number of polynomials in trim problem occurs in the displacement and 
mixed methods, as will be discussed later. Table 1 summaries the number 
of harmonics and required number of polynomials for the trim case. It is 
noteworthy that the unknown trim parameters 6 0 , 0 C and $ t are included in 
the aerodynamic force term F(t) and should be moved to the left hand side 
in the matrix formulation for solution. 



DISPLACEMENT METHOD 


For the displacement formulation, we start with Hamilton’s law of varying 
action 

/"' ( 61 + 6q r Q„ )it - £ |!' = 0 (9) 

Jti i-1 0< li 

With a weak constraint of momentum [19], the Hamilton’s weak principle ( 
primal or displacement form) is then 

6 /“' Ldt + f* SWnc dt - £ 6 gi . Pi I If = 0 (10) 

where the trailing term - £"=i 8q, ‘Pi can be thought of as the virtual action 
6 A entering (or leaving) the time boundaries f, and tj. Thus, one might 
interprete equation (10) as a variational statement of elasto-dynamics, “the 
variation of the action plus the virtual action over the time interval u < o < 
tj plus the net virtual action across the time boundaries must sum to zero”. 

Before proceeding further in displacement development, we integrate equa- 
tion (10) by parts for the linear flapping problem which yields 

r (0 + C0 + Kfi - F)6f3di - (p - 0)60 |2> 0 (11) 

Jo 

Obviously, for equation (11) to be valid for all 8@{F), the following equations 
must hold. 


0 + C0 + K0 

= F 

(12) 

P o 

= m 

(13) 

PT 

= ${T) 

(14) 


The above relations imply that the governing differential equation and bound- 
ary conditions for momenta are satisfied by equation (10). Even for a con- 
servative system, initial value problems are not-self-adjoint and thus they do 


not yield the minimum of a functional in the classical sense. However, the 
variational statement is valid. One can think of Eq. (10) as the time integral 
of the virtual work where the virtual displacements at the boundary times 
( Sqi at U and tj) may or may not be set to zero. 

As shown in Ref. [16], certain geometric boundary conditions must be 
satisfied for the admissable set of trial and test functions to have a stable 
numerical formulation. For the trial functions, geometric boundary condition 
on /3(0) or /3(T), must be exactly satisfied. For the test functions 6/3(0) or 
6/3(T) must be zero if either p 0 or pr is unknown, respectively. 

Eq. (11) can also be interpreted as the sum of the Galerkin weighted 
residuals set equal to zero. For an approximate solution /3, it can also be 
interpreted as an error functional £(/3, 6/3). From the concepts of finite ele- 
ments, one can set E = 0 for some limited class of 6/3. However, the error 
functional is not bilinear if either po or pr is unknown. It follows that if 
either po or pr is unknown, 6/3(0) or 6/3(T) must be set to zero to eliminate 
the unknowns from the trailing terms. Thus, for boundary value problems, 
if /3(0) and pr axe prescribed (p 0 unknown) then we set 6/3(0) = 0. On the 
other hand, for initial value problems, if /3(0) and po are known but pr un- 
known, then we set 6/3(T) = 0 but 6/3(0) ^ 0. If 6/3(0) ^ 0, then setting 
E(/3,6/3) = 0 ensures that po — *■ /3(0) which provides weak convergence to the 
natural boundary condition. Similar arguments holds for other combinations 
of end conditions. 

For the periodic problem, /3(0) = /3(T) is a strong condition whereas 
/3(0) = /3(T) is a weak condition enforced through po = pr which can be 


thought of as a Lagrange multiplier. Equation (10) can then be written as 


[ 2 '{060 - C$6/3 - K060 + F60)dt - \ p [60{2v) - 60(0)] = 0 (15) 

Jo 

Since A p is unknown, we must set 60{ 2ir) = 6/?(0) for a bilinear formula- 
tion. For numerical solution, we choose the integrals of Legendre polynomials 
Pj(x) as shape functions for 0 and 60 

0 = ^2 *j(v)9j » 60 = H *i(v)6n ( 16 ) 

i=i «=i 

where N is the number of terms in 0 or 60 and 

* j+ , = V2^ i = 2,3,...,AT-l 

If the Lagrange multiplier is retained as an unknown (not bilinear), then 
there is no constraint on 5f3 and we have JV + 1 equations in A r + 1 unknowns. 
(Recall that the (N + l) st equation is the required constraint (3(0) = /?( 2tt).) 
However, if 8(3 is constrained, one equation is lost and A p is eliminated yield- 
ing N equations in N unknowns. 

• Momentum in the Displacement Method 

Momentum in the displacement method can be obtained in three ways : 

(i) Derivatives of displacements (or rotations) 

PCh) = Mv) = £ ( 1? ) 

i = 1 * 

»here (') = § = Hg = H )', and , = ^ - 1. 

(ii) Using Lagrange multiplier with derivatives of displacements involved 
in integral formulation 
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P(t x ) = \ P + [ 1 \f - 03 - K(3)di (18) 

Jo 

(iii) Using Lagrange multiplier without derivatives of displacements in- 
volved in integral formulation 

Integration of equation (18) by parts yields 

P(t x ) = A P - 03 ft + /V + C/3 - KP)di (19) 

Jo 

The reason that equation (17) is not advisable, is that there is sensitivity 
problem in convergence of derivatives of shape functions, which may 
occur in (i). By virture of the Lagrange multiplier, strategy (ii) (which 
was used in Ref. [13]) improves the accuracy; and strategy (iii), applied in 
this work, gives the best results. As the number of elements and order of 
polynomials increased, (ii) and (iii) approach the same results rapidly. 

With multiple elements M having (N - l) th order polynomials, the flap- 
ping angle (3 is continuous between elements. This requires M — 1 constraint 
equations, and the periodicity provides one more. This is balanced either by 
M unknowns (which are the momenta at nodes) or by the M constraints on 
6/3 required to eliminate those momenta. These nodal momenta can be ex- 
tracted by discretization of the trailing terms of the displacement formulation 
(15). 

Trailing terms = Am[5/3* 1 ^(0) — 2t }] 

+ ( 20 ) 

m=l 

Thus, the A m can be retained or the 6/3 be restricted. The trim constraint 

equations for multi-elements with ( N — l) 1 * order polynomial are 

M -i N 

S / £ cos< ~ dr] = 0 

m=l J ~ l j - 1 
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where 


M -i N 

£ / sin t dr] = 0 

m=l J ~ 1 j-\ 


£ 1 4 m) £ i<| + | *> f)*; + 2 J, m’ rin 2f *,] <*, 

m=l j= 2 

— 4Cro 


0*jr ^ 

~«i + 5^*) „ 


-7T — Ax 


( 21 ) 


t 


i m + ^Af 

2ir, .» , 1 +t? 2x 
1 + 2 M 


(-!<*/<!) 


Numerical difficulties arise in trim problems if too few terms are taken. 
With a one element model, the number of polynomials has to be at least four 
since we have three trim constraints and one periodicity constraint. With 
multiple element modeling, additional periodicity constraints are required 
between elements, thus the need is for more than four polynomials. The 
relations between number of polynomials and required number of elements 
are shown in Table 2. 


MTXED METHOD 

The momenta, like displacements, are treated as additional field variables 
in the mixed method. The mixed formulation derived from Hamilton’s weak 
principle [5] is 

f * (&q T P - bp 7 q - 6H + Sq T Q)dt = (6q T p - 6p T q) (22) 

Jti 

where q denotes the set of generalized coordinates, Q the set of nonconser- 
vative generalized forces applied to the system, L(q,q,t) the Lagrangian of 
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the system, p = dL/dq the set of generalized momenta, and H = p T q — L 
the Hamiltonian. 

For the flapping problem, 



Q = F - CP 


B - 5 * + \ Kf 

P = P 

The mixed formulation for the flapping problem is then 

f (SP p — Sp P — p Sp — KP 6P + F Sp — C p 8P) dt 
Jo 

= (6Pp-6pP)\l* (23) 

where the variations of P and p are independent, and all boundary conditions 
are of the natural (or weak) type. Since the derivatives exist only in virtual 
displacements SP and virtual momentum 6p, it is then possible to imple- 
ment C° continuous shape functions for the test functions and discontinuous 
functions for the tried functions. 

The C° shape functions are chosen to be hierarchical “bubble” functions 

N - 2 

Sp = <5pi(l - r) + 6p 2 T + (1 - t)t ^2 (o < T < 1) 

t=l 

N - 2 

sp = - r) + Sq 2 T + (1 - r)r ^ Sq*Bi(r)fi, (0 < r < 1) (24) 

»=i 

where Bi(r ) = G,_i(P, Q\t) are Jacobi polynomials. P , Q are parameters 
and fi the normalizing factor; and N is the number of terms in Sp or Sp. 

The approximate values of p and P are taken as continuous function within 
the element while allowing for distinct, discrete values on the boundaries of 



14 


elements [15]. 

V = £ Pj Q j( r )> £ = £ & Q i( r ) (0 < r < 1) 

j=i 

P = Pi, 0 = Pu (r = 0) (25) 

P = P2, P = @2, ( r = 1) 

Note that there are JV — 1 p terms and one unknown p giving N unknown p 
(or (3) terms. The 6/3 and 6 p are assumed as free variations which exist even 
when p and (3 are prescribed. Hodges and Hou [9] showed the derivation of 
shape functions aj(r) and normalizing factors fi for mixed p- version finite 
element. 

Unlike in the displacement method, the trial functions /?, p are discontin- 
uous between elements in the mixed method, while the test functions 6/3, 6 p 
are continuous between elements including the periodic connection. The rea- 
son is to be found below. 

First, we rewrite the trailing terms of the mixed formulation (21) for a 
one-element, model 

(60 p — 6p 0) 1 2 0 '= \ p [60(2 *) - 5/3(0)] - A^[5p(2x) - 5p(0)] (26) 

where 

Pi = h = A p , 0\ — 02 — Ap 

Obviously, the continuity of 60 and 6 p is required to eliminate the unknowns 
A p , A 0 at end of the period. 

For hp finite elements, the nodal momenta can be recovered by follow- 
ing the same steps as equation (20) and by discretizing the time domain in 
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trailing terms 

Trailingterm = 6/3^(0)]+ ^^ m+1 ^m+i)] 

m=l 

M- 1 

-A«Vp ,M, (2-r)-^ m (0)]- ■£ Aj">[« P ( ” , (t„ +1 )-«p"" +,) (t„ + .)] (27) 

m=l 

Thus, the nodal values Aj, m) , A^ m) can then be found by solving equations 
associated with independent coefficients (i = 1,2, . . . ,N — 1). 


The trim equations for the mixed method are 



where 

t = t m + tA< 

Again, there is a minimum requirement, for the number of polynomials in 
the mixed method trim problem. With Lagrange multiplier, the periodicity 
constraints on generalized coordinates are enforced “‘implicitly”; thus, the 
minimum polynomial number for trim is six for any number of elements 
since there are two generalized coordinates, displacement and momentum. 


This relation is depicted in Table 3. 
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• Flap-Lag 

We use the same model as in Ref. [8] for the flap-lag problem. The 
schematical model is shown in Fig. 1. It consists of a centrally hinged rigid 
blade, with flap and lag restraint springs to simulate the elastic coupling 
characteristics of the actual blade. The spring stiffnesses are chosen so that 
the uncoupled, rotating flap and lag natural frequencies coincide with the 
corresponding first mode rotating natural frequencies of the elastic blade. 
The effects of induced flow, blade equilibrium, and aerodynamic stiffness on 
aerodynamic terms are retained. We further assume that the drag coefficient 
Cdo is small with respect to the lift curve slope, and reversed flow is neglected. 
The nonlinear equations of motion of this model axe [7] 

j3 + sin/3cos/3(l + C) 2 + (P — !)(/? — Ppc) + = J o Fpfdf (29) 

cos 2 /3C -2sm/3cosP{l + C)0 + WC + Z{0-0 pc ) = cos /3 F ( r dr 

where Fp and are the nondimensional airloads perpendicular to the blade 
and respectively in the flap and lead directions. The elastic terms approxi- 
mated by small (3, ( formulas [12] and small pitch angle assumption are 

P = l+u;l=p 2 , W = u>l Z = R{u\-J 1 (i )e (30) 

For the linearized model, we make small perturbations about a periodic 
equilibrium motion of the nonlinear system. Denoting 0(f), C(f) /9(f) 
as the equilibrium values of 0, £ and (3 and one can write the perturbation 
expansions : 

C = C + *C, (3 = P + 8/3, O = e + 9 0 6(3 + e c 6( (31) 

where 

0 = 0 O + 0, sin f + 0 C cos t + 0p(j3 — (3^) + 0<C 
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and Op and 0 ( are pitch-flap and pitch-lag coupling parameters. Positive Op or 
8( implies that the blade pitches up due to positive flapping or lead inplane 
motions. 

Substituting Eq.(31) into the nonlinear equations (29), cancelling the 

appropriate equilibrium terms and neglecting the higher order quantities 
- - - - - 2 

C, C, P, 00, and 0 with respect to unity; we obtain the linearized equation 
for flap-lag forced response : 


{3} + p®i{8) + ™i{k}- 

[£<<)] { \ } + (Af(f)] | l } + [W(()| | c $ /a J + {0(f)} ft, (32) 

where 

1 + + 0 ) 

— 2j0(l + |/isf) + 20 

-2^+| Iict0 + i$) \Q{4> + \ pet 0 + 0) — 200 

+ ^( 1 + 3 Mt) - 20 +2 a i(l + l^ <_ ) 


[C(t)) = 


(33) 



p 

+|(f ticE+2p*sid) _ 
— 8^(1 2^i 2 5 2 F) 


Z -2$ - 2ifici(l<f> + §/ 3 ) 
+^(|/xcf + 2^sict) 
—lj$fi*c 2 i + — /3pc) 

■(wj - wj) 

+^(|/icf + 2// 2 5tcf)] 


Z + ^ C f(f*+f£) 
-2^(|/icf+2/i a «fcE) _ 
+2 |/3(/i 2 c 2 i- n 2 s 2 t- Ifist) 

- 19 C ( 1 + f/isf + 2/iVi") 

+/jct0 

-2fi 2 cH) 

+2/MFj8(f$+|£ + 2/ie#)] 

+ J*cl$(l + |/X5t) 

H-^( + 2 fi 2 stct) 

+£(i + |m 5 0] 

( 34 ) 
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[L(i)} = 


[M(i)) = 

[N(t)} = 


W)} = 


-1(1 + | si) -P - J(f^cT+ ifsid) ■ 

l(24> + I fact + £) *(3&ict + 2 fcHp) 

-20(1 + | fist) + 20 + 2fi 2 ctsi) 

-Z - R(u>> - *1)00 

’ 2fi 2 s 2 t) -J(l + f/«*)" 

. -**! + §/»-*) i(f*). . 

-1 0 

.0 -|(1 + |/xsf+ 2/x 2 5 2 <) 

I 1 

k Z + R{^1 — J 


(35) 


(36) 


(37) 


(38) 


Due to the nonlinear effects of flap-lag coupling for trimmed flight with 
constant Ct , iterations for the solutions are required. To speed up the con- 
vergence, formulas for 9o,0, and 0 C in Ref. [ 11 ] by harmonic balance can be 
used as the first guess. 


ERROR ANALYSIS 

The shape functions are integrals of Legendre polynomials, $( 77 ) for the 
displacement method; and are Jacobi polynomials Bj(r) for the mixed method 
These polynomials can be related as follows 

= Cj (»? 2 - 1) Bj. 1 , j> 2 (-1 < v < 1) (39) 

where r = a ^- and Cj is constant. 

In order to express the characteristics of shape functions; errors are mea- 
sured at the zeros (roots) of shape functions. The roots of integral of Legen- 
dre polynomials, for example, are more closely grouped near the endpoints. 
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This indicates that the polynomials are able to oscillate with increased fre- 
quency near the endpoints, and are better suited for approximating singular 
behavior which occurs at the endpoints than are trigonometric functions or 
uniform /? -version meshes. For trigonometric functions, the roots are evenly 
distributed in the interval -1 < x < 1 and they have the important property 
that they “absorb” singularties at the finite element boundaries. Equation 
(39) also indicates that the zeros of Jacobi polynomials and the integrals of 
Legendre polynomials are equivalent. Table 4 and 5 show the number of 
zeros and polynomial order in the mixed and displacement formulations. 

The shape functions in the displacement method start with linear func- 
tions ( N = 2), since the order of the polynomial n = N — 1. We thus need 
2 points in each element for error analysis, which would be the nodes of el- 
ements. For IV = 3, the quadratic polynomial needs an additional point to 
measure error which is the zero of next higher order polynomial. By keeping 
the state variables in their primative forms, the mixed method is able to start 
with constant (TV = 2 in Eq. (25) ) for its shape function since the order of 
polynomial n — N — 2. Following the same argument as in the displacement 
method, the measured points for N th polynomial are the zeros of (N + l) tJl 
order polynomial. 

Denoted by the norm error for displacement and mixed methods 

are defined seperatly as 


E d = 


1 


e»=i {(A/ji">)’ + + [(ApM)’ + ES^Api"’) 2 ] • (WY) 

E" 1 {(/S' 1 ) 1 + E"?^” 1 ) 2 + [(pi? 1 ) 1 + i£V(pS~V1 • W 2 } 


e„ = 


Em-1 {(A^"))= + £? ,(A^"V + [(Ap(”>)* + E* ,(Ap} m) ) 2 ] ■ (W)*} 

i E" , {(/S 0 ) 2 + E Uifc'r + [(pS 0 ) 2 + E" 2 (p£’) 2 J • (WO 2 } 


where weight W = T/ 2tt; A0 {m) is the error at node of m th element, and 
Afij is the error at j th zero within element. 

For linear flap quasi-trim, trim and linearized flap-lag quasi-trim prob- 
lems, essentially exact solutions can be obtained by Fourier series analysis. 
Since the Fourier series approach is a special case of time finite elements 
using one element, Figures 2, 3 and 4 illustrate the convergence behaviors 
of three approaches for quasi-trim and trim cases by single element model- 
ing. Fourier analysis obviously gives excellent accuracy with respect to the 
number of floating point operations and CPU. Realizing that even the mixed 
formulation has almost twice the number of unknowns as in the displacement 
formulation and the Fourier analysis, and thus increase the number of oper- 
ations and the size of the matrix; it yields a very sparse matrix which would 
be significantly influence the computation efficiency as will be shown later. 

With multi-elements modeling, Fig. 5 (in which number of elements 
progresses as 1,2,4,8,16 and 32) compares the number of floating-point oper- 
ations (see formulation in Ref. [13]) between the displacement method and 
the Fourier analysis. It shows that the displacement method is not compet- 
itive with classical Fourier analysis. Figures 6 through 9 compare the norm 
error with regard to the number of degrees of freedom for quasi-trim and 
trim cases. Without consideration of computational efficiency, the Fourier 
approach always gives the best results. However, the main objective here 
is to obtain satisfying accuracy with highest efficiency (or minimum CPU). 
In Figures 10-13, norm errors are plotted against CPU (on UNIX Sequent 
S81 computer) for the linear flap problem; and Figures 14 - 17 are for lin- 
earized flap-lag problem. Figures 10, 11 show that the mixed method is very 
competitive with the Fourier approach and is faster than the displacement 



method for a specified error requirement for the quasi-trim case. The advan- 
tage of mixed formulation in computation efficiency becomes obvious for the 
trim case as shown in Figures 12, 13. It is shown that at any given error 
criteria or CPU, the mixed method provides an optimum choice of number 
of elements and polynomial order (i.e., hp trade-off) due to the advantage of 
matrix sparsity. 

In general, 4 linear elements in the mixed method is appopriate for mod- 
erate error criteria (1%), and 4 to 8 higher order elements for more stringent 
error criteria. However, the hp trade-off is not obvious in the displacement 
method, and the family of curves tends to converge with a slope of -5 to 
-6. More precisely, the relation between the norm error and CPU in the 
displacement method can be formulated as following 


e 2 .CPU, ^ 

% - { cpu; ) 


(40) 


In another words, for a factor of 10 improvement in accuracy, displace- 
ment method needs 1.5 times the CPU; and a factor of 100 improvement 
needs 2.5 times the CPU. Comparing both methods for the flap-lag problem 
in Figures 14 - 17, we can see clearly that mixed method is much faster than 
the displacement method for any given error criteria. Furthermore, the spar- 
sity of matrix for mixed method has not yet been fully utilized, which will 
increase the advantage of the mixed method. 


CONCLUSIONS 


Finite elements in time (FET) are applied on rotorcraft trim problems 
in which Fourier series analysis, displacement and mixed methods are im- 
plemented and verified for accuracy and efficiency. Examples include linear 



flap and linearized flap-lag problems. It is proved that Fourier series can be 
thought of as a special case of finite elements in time, and the results also 
show that the mixed method in FET can outperform conventional Fourier 
analysis for rotor trim problems. The shape functions for different methods 
are chosen according to the characteristics of associated formulations. Thus, 
the integrals of Legendre polynomials are preferrable for the displacement 
method while Jacobi ploynomials are preferrable for the mixed method. 

The mixed method shows the advantage of efficiency at a given accuracy. 
From hp trade-off analysis, it seems that 4 to 8 elements in the mixed method 
provide best results for any order polynomial for any given error criteria, 
while the hp trade-off is not sensitive in the displacement method, for which 
all the curves have similar convergence behavior. 
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Table 1. No. of Polynomials for Trim Using Fourier Senes 


No. of 


Harmonics Nh 


Fourier Series 


0 No. of Poly | Req. Poly No. 

= 2N h + 1 


singular matrix 


Table 2. No. of Polynomials for Trim Using Displacement Method 

Displaceme nt Method 

M N 1 0 I No. of Poly I Req. Poly No. 


n 


— X 

— 


= MN 
2 * 

3* 

4 

4* 

6 

6 

8 


= M + 3 
4 
4 

4 

5 

5 
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* singular matrix 
















Table 3. No. of Polynomials for Trim Using Mixed Method 


^ it 

Mixed Method 

M 

N 

0 

P 

No. of Poly 
= 2M(N — 1) 

Req. Poly No. 

1 

2 



2* 

6 

1 

3 



4* 

6 

1 

4 



6 

6 

2 

2 



4* 

6 

2 

3 



8 

6 

3 

2 



6 

6 

* singular matrix 


Table 4. No. of Zeros and Polynomial Order of Displacement Formulation 


N 

Description 

No. of Zeros including 
end points 

No. of Zeros 
in intervals 

2 

Linear 

2 

0 

3 

Quadratic 

3 

1 

4 

Cubic 

4 1 

2 

5 

Quatric 

5 

3 


Table 


N 

Description 

No. of Zeros including 
end points 

No. of Zeros 
in intervals 

2 

Constant 

0 

0 

3 

Linear 

1 

1 

4 

Quadratic 

2 

2 

5 

Cubic 

3" 

3 

6 

Quatric 

4 

4 
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Fig. 1 Centrally hinged, spring restrained, rigid blade representation. 
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Fig. 2 Norm error vs. No. of operations for flapping (quasi-trim) 
with one element model. 
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Fig. 7 Norm error vs. degrees of freedom for flapping (quasi-trim) 
by mixed method compared with Fourier approach. 
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Fig. 8 Norm error vs. degrees of freedom for flapping (trim) 

by displacement method compared with Fourier approach. 
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Fig. 9 Norm error vs. degrees of freedom for flapping (trim) 
by mixed method compared with Fourier approach. 
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Fig. 10 Norm error vs. CPU for flapping (quasi-trim) by 

displacement method compared with Fourier approach. 
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Fig. 11 Norm error vs. CPU for flapping (quasi-trim) by 
mixed method cpmpared with Fourier approach. 
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Fig. 12 Norm error vs. CPU for 
ment method compared w 
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Fig. 13 Norm error vs. CPU for flapping (trim) by mixed 
method compared with Fourier approach. 
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Fig. 14 Norm error vs. CPU for flap-lag (quasi-trim) by 
displacement method. 
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Fig. 15 Norm error vs. CPU for flap-lag (quasi-trim) by 
mixed method. 
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Fig 16. Norm error vs. CPU for flap-lag (trim) by displace- 
ment method. 
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Fig. 17 Norm error vs. CPU for flap-lag (trim) by mixed 
method. 



